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ABSTRACT 

We study the final architecture of planetary systems that evolve under the combined effects of 
planet-planet and planetesimal scattering. Using N-body simulations we investigate the dynamics of 
marginally unstable systems of gas and ice giants both in isolation and when the planets form interior 
to a planetesimal belt. The unstable isolated systems evolve under planet-planet scattering to yield 
an eccentricity distribution that matches that observed for extrasolar planets. When planetesimals 
are included the outcome depends upon the total mass of the planets. For M tot > 1 Mj the final 
eccentricity distribution remains broad, whereas for M tot < 1 Mj a combination of divergent orbital 
evolution and recircularization of scattered planets results in a preponderance of nearly circular final 
orbits. We also study the fate of marginally stable multiple planet systems in the presence of planetes- 
imal disks, and find that for high planet masses the majority of such systems evolve into resonance. A 
significant fraction lead to resonant chains that are planetary analogs of Jupiter's Galilean satellites. 
We predict that a transition from eccentric to near-circular orbits will be observed once extrasolar 
planet surveys detect sub- Jovian mass planets at orbital radii of a ~ 5 — 10 AU. 
Subject headings: solar system: formation — planetary systems: protoplanetary disks — planetary 
systems: formation — celestial mechanics 



1. INTRODUCTION 

Different dynamical mechanisms are commonly in- 
voked to explain the architecture of the outer Solar 
System and extrasolar planetary systems. In the So- 
lar System, scatte ring of small bodies ( " planetesimals" ) 
by the ice giantsjFernandez fc Id 11984 : llda et al.ll2000i 
iKirsh et al.ll2009f ) is thought to drive outward planetary 
migration and concomitant capture of Pluto and other 
Kuiper Belt Objects into r esonance (jMalhotral HMl 
iMurrav-Clav fc Chiang] 120051 ) . The effects of planetes- 
imal scattering on the gas giants are small er but still sig- 
nificant, for example in th e "Nice model" (Tsigani s et al.l 
120051 iGomes et all 120051 ) where a divergent resonance 
crossing between Jupiter and Saturn triggers the Late 
Heavy Bombardment. The presence of small bodies 
around other stars can be infe rred from observations of 
debris disks (|e.g. Wvattl I2008T ). but as yet there is no 
evidence for a dynamical role of planetesimals in known 
extrasolar planetary systems. At radii where tidal ef- 
fects are negligible (roughly a > 0.1 AU) the eccentricity 
distribution of extrasolar planets matches relatively sim- 
ple models of gravitational scattering among a system 
of two or more massive planets that typically include 
neither planetesimals nor resi dual gas (jChatteriee et al.l 
2008; Uuric fc Treriiam^[200l . 

The success of pure planet-planet scattering mod- 
els does not imply that other dynamical processes can 
be ignored. The observed distribution of semi-major 
axes of extrasolar planets at small orbital radii re- 
quires the existence of an additional dissipative process 
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( Ada ms fe Lau ghlin 20Q3|), most probably gas disk mi- 
gration (|Lin fc Papaloizoul I1986D. which will itself af - 
fect planetary eccentricity (M oorhead fc Ada ms 2005). 
At larger orbital radii simple arguments suggest that 
a dynamically significant external reservoir of planetes- 
imals ought to be a common feature of young plane- 
tary systems. The formation of gi ant planets becomes 
increasingly difficult at large radii ( Poll ack et al.l [19961 : 
iKokubo fc Idall2002f ). and hence it is probable that disks 
of leftover debris surround the zone of giant planet forma- 
tion in most young systems. The typical masses of plan- 
etesimal disks are unknown, but values of 30-50 M® that 
are comparable to those inferred for the early outer So- 
lar System are plausibly typical, since they are consistent 
with disk masses estimated from astronomical observa - 
tions of the youngest stars ([Andrews fc Williams! [2005) . 
The dynamical effect of such disks on currently observed 
extrasolar planetary systems would be small, since radial 
velocity surveys preferentially detect planets that are ei- 
ther massive (and hence largely immune to influence from 
planetesimal disks) or orbit at very small radii where the 
mass of leftover debris is negligible. 

Pooling knowledge from the Solar System and extraso- 
lar planetary systems motivates consideration of a model 
in which planet formation typically yields a marginally 
unstable system of massive planets in dynamical contact 
with both a residual gas disk and an exterior planetesi- 
mal disk. In this Letter we ignore the gas disk and study 
the subsequent evolution under the combined action of 
planet-planet and planetesimal scattering. We do not 
attempt to model the full distribution of extrasolar plan- 
etary properties (which would require the inclusion of 
hydrodynamic effects), but rather focus on how plan- 
etesimal disks affect the final eccentricity of extrasolar 
planets at moderately large orbital radii. 

2. METHODS 



2 



Raymond, Armitage & Gorelick 



We assume that the gas-dominated epoch of planet 
formation is sufficiently distinct from the subsequent 
phase of planet-planet and planetesimal scattering that 
it makes sense to study the latter with pure N-body sim- 
ulations. We focus on two large ensembles of runs. The 
highmass set comprises 1000 integrations of three planet 
systems in which the masses of the planets are chosen 
randomly in the range Mg at < M p < 3Mj, with a distri- 
bution, 



which matches that observed (jMarcv et al.l [2008V The 
observed distribution is derived from an incomplete sam- 
ple that represents (in the context of our model) the 
distribution after scattering, but these subtleties do not 
matter for our purposes. The lowmass set is identical ex- 
cept that we sample a wider swath of the mass function 
between 10M® and 3Mj. The planets are initially placed 
in a marginally unstable configuration defined by circu- 
lar, nearly coplanar orbits with a separation of 4-5 r^. m , 
where the mutual Hill radius. 
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Here ai and ai are the planets' semi-major axes, M\ 
and M2 their masses, and M* is the stellar mass. With 
this spacing the instability t i mescal e is relatively long 
([Chambers. Wetherill fc Boss! 119961: iChatteriee et al.l 
I2008h (the median timescale before the first planet-planet 
encounter was 0.3 Myr for the highmass integrations 
without disks). Our initial conditions are only a small 
subset of the architectu res predicted from giant planet 
formation models (|Thommes. M atsumu ra fc Rasiol 
I200a iMordasini. Alibert fc Benzll2009D . though broadly 
consistent with scenarios in which giant planets are 
captured into mean-motion reso nances during t h e late 
stages of gas disk evolution (jThommes et al.1 I2008T ) 
prior to being removed from resonance by t urbulent 
perturbations (|Adams. Laughlin fc Bl och 2008). Each 
integration is repeated twice, once with just the three 
planets 5 and once with an external planetesimal disk 
whose inner radius of Oi n = 10 AU is 2 Hill radii beyond 
the orbit of the outermost planet 6 . The inner edge of 
the disk lies within the radius where a test particle 
in the restricted 3-body problem would be stable, so 
the disk is in immediate dynamical contact with the 
outer planet. The planetesimal disk is represented by 
1000 bodies distributed between 10 and 20 AU with a 
^disk oc r _1 surface density profile and a total mass of 
50M®. 

We integrate t h ese s ystems using the MERCURY 
code ([Chambers! 1 19991) for 100 Myr. The intc- 
rator uses the symplect ic Wisdom-Holman mapping 
Wisdom fc Hol man 199j]) for well-separated bodies, and 
the Bulirsch-Stoer method when objects are within N 

5 Note that the simulations without plane tesimal disks are iden- 
tical to the mixedl and mixed2 cases from (Raymond et all 120081 
[200911 . 

6 In the Nice model, a separation of ~ 3 — 4 Hill radii between 
Neptune and the outer planetesimal dis k is needed to mat ch the 
timing of the Late Heavy Bombardment ( Gomes et al. 2005j). The 
spacing of 2 m means that our models evolve on a somewhat 
shorter time scale. 
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Fig. 1. — Cumulative eccentricity distributions for observed 
extra-solar planets (thick grey line; lighter grey for minimum 
masses M < Mj) as compared with distributions from our sim- 
ulations. The model distributions include only the eccentricity of 
the innermost (and hence most easily detected) planet. 

mutual Hill radii, where N = 3 for our case. Planets 
were removed if their orbital distances were smaller than 
0.1 AU ("hit Sun") or exceeded 100 AU ("ejection"). 
Collisions were treated as inelastic mergers conserving 
linear momentum. 

A large ensemble of simulations includes some cases 
that are much harder to integrate accurately than the 
majority. To make the best use of our computational 
resources we adopted a default timestep (20 days) that 
results in accurate integrations (as measured by the frac- 
tional orbital energy conservation dE/E) for the typical 
case. We then identified those runs (about 10%) in which 
energy was not adequately conserved and re-ran them 
with a smaller timestep. For runs without disks we re- 
ran cases with dE/E > 10~ 4 with a timestep of 5 days, 
while for the runs with disks we re-computed cases with 
dE/E > 5 x 10 -4 with a timestep of 10 days. A small 
number of the re-run simulations (typically 15-35) still 
did not meet our energy criterion and were discarded. 

3. RESULTS FOR MARGINALLY UNSTABLE PLANETARY 
SYSTEMS 

In the absence of planetesimal disks our model plane- 
tary systems are typically unstable on Myr time scales. 
There are also systems that are stable over the 100 Myr 
duration of our runs. In our initial analysis we assume 
that the typical outcome of giant planet formation is a 
system that, in the absence of a disk, would be unsta- 
ble. We therefore analyze the subset of disk-less simula- 
tions that are unstable, and compare the results to the 
matched set of simulations that include disks. This is not 
a perfect one-to-one comparison, since the chaotic nature 
of the evolution means that disk-less planetary systems 
can display different instability time scales in the pres- 
ence of even negligible perturbations. Nonetheless we 
do observe statistical differences between the evolution 
of systems (with disks) that correlate with the stability 
of the disk-free systems, and hence it makes sense as a 
first approximation to separately consider the results for 
stable and unstable cases. 

The results of our disk-less simulations agree with 
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prior studies (iChatteriee et alj l2008t Uuric fe Tremainel 
120081 iFord fc Rasiol I2008D . Scattering from initially 
unstable initial conditions frequently leads to the loss 
of one or more planets via ejection or collisions and 
sets up a broad eccentricity distributi on (Rasio fc Ford 
ll996HWeidenschilling fc Marzaril[l99l iLin fc Idalll997l ). 
Scattering among equal-mass planets produces larger 
eccentricities than scattering of unequal-ma ss planets 
(jFord. Rasio fc Yul[200a iRavmond et al.ll2008fl . Figure Q] 
compares the final eccentricities obtained from the unsta- 
ble highmass simu l ations and the o bserved distribution 
([Butler et all 120061 : ISchneider![2009] ) of extra-solar plan- 
ets. They are in good quantitative agreement. The ec- 
centricity distribution from the unstable lowmass sim- 
ulations without disks is shifted toward lower values 
(|Ravmond et all 120081 : [Ford fc Rasiol 120081 ). and fits the 
obse rved distribution of extra-solar planets with M p < 
Mj ([Wright et al.l 12008( 1. Our model therefore exhibits 
evolution that is consistent with current observations of 
extrasolar planetary systems, which as we noted previ- 
ously are mostly of systems at such small radii that plan- 
etesimal disks are dynamically unimportant. 

At larger radii we expect that both disks and planet- 
planet scattering will play a dynamical role. A wide 
range of outcomes is then possible. Exchange of 
energy and angular momentum between the planets 
and the planetesimal disk leads to planetary migra- 
tion (iFernandez fc Ip||1984tlMurrav et al.lll998tllda et all 
120001 iKirsh et al.l l2009l )7which can be either stabiliz- 
ing or destabilizing. A low mass planet adjacent to 
the disk scatters planetesimals inward, resulting in di- 
vergent migration that is often stabilizing unless reso- 
nance crossing excites eccentricity to the point of trig- 
gering instability. Alternatively, an outer massive planet 
interacting with the disk directly ejects planetesimals 
and migrates inward, compressing the system and lead- 
ing to instability or resonant capture. An equally im- 
portant effect is that the disk can act to recircularize 
the orbits of scattered planets after dynamical instabili- 
ties (Thommes. Duncan fc Levisonll 19991 : [Ford fc Chiang! 
12007 ). To illustrate how significant recircularization can 
be we ran a small additional set of idealized experiments 
in which a single planet with mass M p on an orbit with 
a = 10 AU and e — 0.5 begins to interact with our initial 
planetesimal disk. For 10 4 — 10 5 years (longer for smaller 
M p ), e is damped roughly exponentially with a damping 



time scale t e , defined via, 
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of t e w 5 x 10 4 years, independent of M p . The sub- 
sequent evolution was highly mass-dependent: for low- 
mass planets, e continued to decrease on much longer 
timescales (0.36, 0.63, and 4.6 Myr to reach e < 0.1 
for M p = 10M e ,30M ffi , and Ms, respectively). Massive 
planets disrupted the disk, halting dynamical friction. 
The total decrease in e for M p > Mj was 0.15 or less, 
corresponding to an increase of 1.5 AU or less in perihe- 
lion distance. 

Figure [2] illustrates the diversity of outcomes from our 
simulations that include planetesimal disks. We split our 
simulations into three mass bins (the Solar System's gi- 
ant planets fall into the middle bin) and three stability 



categories. In "stable" systems there are no close encoun- 
ters between planets and no large-scale change in system 
architecture (the ordering of the planets is preserved and 
all planets survive). "Moderately stable" systems experi- 
ence substantial perturbations - which may be due to res- 
onance crossing in high-mass systems or close encounters 
in low-mass systems - that are nonetheless insufficient to 
alter the architecture. "Unstable" systems undergo close 
encounters leading to architectural change. 

Subsets of our runs show dynamics analogous to that 
studied for the Solar System and for extra-solar plan- 
etary systems. At high masses planetesimal disks sta- 
bilize about 30% of cases but planet-planet scattering 
leading to the loss of one or more planets is still com- 
mon. Quantitatively, the median eccentricity is reduced 
(Figure [l} but many highly eccentric systems remain. 
As planet masses decreases the dynamical importance of 
planetesimals grows. For low-mass systems, the masses 
of the planets and the planetesimal disk are comparable 
and planetesimal scattering inevitably leads to migra- 
tion. Divergent crossing of mean-motion resonances (one 
example of which is shown in the center panel of Figure^ 
can result in abrupt changes to planetary semi-major axis 
and eccentricity that qualitatively resemble those seen i n 
the Nice model ([Tsiganis et al.ll2005tlGomes et al.ll2005D . 
We also see behavior that resembles an alternative So- 
lar System model in which Uranus and Neptune formed 
in the Jupiter-Saturn region and wer e scattered outward 
(jThommes. Duncan fc Levisonll 1999t ) (top left panel). At 
the lowest masses even highly unstable systems rarely de- 
stroy any planets because recircularization of scattered 
planets is efficient (top center panel), though re-ordering 
of planets is common. In summary, dynamics character- 
istic of the outer Solar System is common among low- to 
medium-mass planetary systems. 

The main prediction of our model is the statistical dis- 
tribution of planetary eccentricity as a function of planet 
mass. Figure [31 shows the distribution of the eccentricity 
of the innermost surviving planet as a function of the to- 
tal mass in surviving planets. In the absence of disks we 
observe similar behavior across all system masses - the 
shift to smaller eccentricities for the lowmass runs, seen 
in Figure [TJ is not visually apparent. When disks are 
included, the eccentricity distribution divides into two 
distinct regimes: a low mass regime in which planetesi- 
mal dynamics dominates to yield low eccentricities and 
a high mass planet scattering dominated regime where 
planetesimals play a minor role. For our specific param- 
eters (inner edge disk edge at 10 AU, and a disk mass 
of 50 M s assumed to be typical for a stellar metallic- 
ity Z = Zq) the transition between these regimes occurs 
for system masses M cr - lt ~ 1 Mj. We predict that sys- 
tems whose giant planets orbit between 5 and 10 AU, and 
which have a total mass below 1 Mj, will typically have 
low eccentricity orbits. This critical mass should scale 
roughly linearly with the stellar metallicity, as we expect 
the initial planetesimal disk mass to be proportional to 
Z . We expect the same qualitative behavior even if the 
zone of giant planet formation e xtends to larger radii 
(|Goldreich. Lithwick fc S ari 2004), though in this case 
the low eccentricity regime would only be observable fur- 
ther out. 
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Fig. 2. — Evolution of a range of planetary systems interacting with planetesimal disks. Each panel shows the evolution of the semimajor 
axis a, perihelion and aphelion distances q and Q for the three planets of a given simulation (planetesimal particles are not shown). All 
simulations have the same x axis scale, but different simulations have different y axis scales, so the region from 5-10 AU is shaded in 
grey for each case. Each column shows simulations in a given mass range: low-mass (total initial planet mass Mtot < 0.5Mj), medium- 
mass, (0.5Mj < Mtot < 2Mj), and high-mass (Mtot > 2Mj). Each row groups simulations by outcome: unstable cases underwent close 
encounters between planets, moderately stable systems underwent significant orbital changes during the simulation but the system remained 
stable, and stable systems did not undergo any close encounters between planets and remained stab le throughout. The evolution in the 
center and top-left panels are qualitativ ely similar to different models of early Solar System dynamics (Thommcs, Duncan & Lcvison 1999; 
ITsiganis et~aL1l2005l : IGomes et aTll2005l) . 

4. RESULTS FOR MARGINALLY STABLE PLANETARY 
SYSTEMS 

Although the eccentricity distribution of extraso- 
lar planets is consistent with the hypothesis that all 
newly formed multiple systems are unstable in the 
absence of disks, this conclusion may also be biased 
by selection effects. Most known extrasolar planets 
probably suffer ed significant gas disk migration prior 
to scattering (iLin, Bodenheimer fe Richardson] |1996|: 
Trillin sTet al.1 119981 : iBodenheimer, Hubickvi fe Lissaueri 




2000), so the high incidence of instability may be a conse- 



0.3 1 
Total Mass (MJ 

Fig. 3. — Final eccentricity of the innermost planet as a function 
of the total mass in surviving planets for the highmass (black) 
and lowmass (grey) simulations. The plotted sample shows those 
systems that were unstable without disks (bottom panel), together 
with the matched sample including disks (upper panel). Disks 
result in a sharp transition to a low-eccentricity regime for system 
masses below about one Jupiter mass. 



quence of migration rather than formation. With this in 
mind we have separately analyzed those (previously ex- 
cluded) systems that were stable in the absence of disks 
to see what impact disks have on them. As expected, 
low eccentricity outcomes predominate. Resonances are 
also common: about 70% of all stable highmass sim- 
ulations and 1/3 of stable lowmass simulations include 
at least one pair of planets in the 3:2 or, more often, 
the 2:1 mean motion resonance. This is a much higher 
probability of resonance capture than occurs for pure 
plane t-planet scattering without disks (|Ravmond et al.l 
2008), and it also exceeds the fraction of resonant sys- 
tems that are expected to survive the gaseous disk phase 
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in the presence of turbulence (lAdams. Laughlin fc Blochl 
2008). Most surprisingly within the highmass set a 
substantial fraction (about 1/3) of stable systems be- 
come locked into mean-motion resonances that involve 
all three of the planets - analogs of the Laplace reso- 
nance among Jupiter's Galilean satellites. Chains of res- 
onances 7 arise preferentially in higher-mass systems and 
in systems where planetesimal-driven migration causes 
compression rather than divergent migration. Detection 
of high mass planets at the relevant radii (between 5- 
10 AU) should soon be possible via astrometric or di- 
rect imaging techniques, and observation of resonant 
chains would be consistent with our model. Intrigu- 
ingly, the recently-discovered triple planet system HR 
8799 may be in a 4:2:1 resonant c hain (|Marois et alj2008t 
iFabrvckv fc Murrav-Clavl [20091 . Determining whether 
capture into resonance was initiated by a gas or plan- 
etesimal disk may be possible via detailed compari- 
son o f the outcome of resonant ca pture in planetes- 
imal (iMurray-Clay fc Chiand l200l versus gas d isks 
(|Lee fc Pealdl2002tlAdams. Laughlin fc Bloch|[200l . 

5. CONCLUSIONS 

Circumstantial evidence suggests tha t many observed 
properties of the outer Solar Syste m (IMalhotr al 11995 ) 
and of extrasolar planetary systems dRasio fc Fordlll996 : 
iWeidenschilling fc Marzarl Il996t iLin fc Idal 119971) may 
be attributable to the dynamical effects of planetesimal 
scattering and planet-planet scattering. Here, we have 
studied the predicted architecture of planetary systems 
that results from the joint action of both mechanisms. 
We have argued that this regime will be relevant once 
lower mass extrasolar planets are discovered at larger 
orbital radii than those currently known. Generically 
we predict that a transition to "Solar-System-like" archi- 
tectures, characterized by near-circular orbits and rela- 
tively stable planetary separations, will be observed once 

7 Our definition of resonance requires one resonant argument to 
librate with an amplitude A < 150° . Roughly half of the resonant 
systems were deep in the resonance, with A < 60°. This fraction 
appears to be independent of the number of planets in resonance, 



surveys detect planets in the regime where planetesimal 
disks play a dynamical role. Our simulations suggest 
that the transition is a surprisingly sharp function of to- 
tal planetary system mass, and that it occurs for system 
masses a factor of several larger than the initial planetes- 
imal disk mass. 

Our initial conditions do not sample anything ap- 
proaching the full range of initial planetary system ar- 
chitectures. We believe that the existence of a tran- 
sition between typically eccentric and near-circular or- 
bits is a general feature of joint models of planet-planet 
and planetesimal scattering, but the transition mass and 
minimum orbital radii at which planetesimal effects be- 
come manifest is of course a function of the poorly known 
masses and radial extent of planetesimal disks. Our re- 
sults suggest that the transition might be seen for sub- 
Jovian mass planets at orbital radii of 5-10 AU, but the 
transition would be pushed to greater orbital radii if giant 
planet formation consumes planetesimals across a wider 
extent of the disk. We also find that the final system 
architecture varies substantially depending on the initial 
separation of the planets. In particular, if planet forma- 
tion yields a mixture of massive systems in initially stable 
orbits, interaction with planetesimals drives a large frac- 
tion of systems into resonance. Whether such systems 
exist should be testable in the near future. 
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as about 1/4 of the highmass resonant chains had A < 60° for both 
pairs of planets. 
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